home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / lisp / kcl / akcl / kcl.lha / c / earith.c < prev    next >
C/C++ Source or Header  |  1987-06-04  |  7KB  |  466 lines
  1. /*
  2. (c) Copyright Taiichi Yuasa and Masami Hagiya, 1984.  All rights reserved.
  3. Copying of this file is authorized to users who have executed the true and
  4. proper "License Agreement for Kyoto Common LISP" with SIGLISP.
  5. */
  6.  
  7. /*
  8.     earith.c
  9.  
  10.     EXTENDED_MUL and EXTENDED_DIV perform 32 bit multiplication and
  11.     division, respectively.
  12.  
  13.     EXTENDED_MUL(D,Q,R,HP,LP) 
  14.     calculates D*Q+R and saves the result into the locations HP and LP.
  15.     D, Q, and R are 32 bit non-negative integers and HP and LP are
  16.     word addresses.  The word at LP will contain the lower 31 (not 32)
  17.     bits of the result and its most significant bit is set 0. The word
  18.     at HP will contain the rest of the result and its MSB is also set 0.
  19.  
  20.     EXTENDED_DIV(D,H,L,QP,RP)
  21.     divides [H:L] by D and saves the quotient and the remainder into
  22.     the locations QP and RP, respectively.  D, H, L are 32 bit non-negative
  23.     integers and QP and RP are word addresses.  Here, [H:L] means the
  24.     64 bit integer (imaginary) represented by H and L as follows.
  25.  
  26.       63 62                  31 30                 0
  27.       |0|0|<lower 31 bits of H>|<lower 31 bits of L>|
  28.  
  29.     Although [H:L] is 64 bits, you can assume that the quotient is always
  30.     represented as 32 bit non-negative integer.
  31. */
  32.  
  33. #include "include.h"
  34.  
  35. #ifdef VAX
  36.  
  37. extended_mul(d, q, r, hp, lp)
  38. int d, q, r;
  39. int *hp, *lp;
  40. {
  41.     asm("    emul    4(ap),8(ap),12(ap),r0");
  42.     asm("    ashq    $1,r0,r0");
  43.     asm("    rotl    $-1,r0,r0");
  44.     asm("    movl    r0,*20(ap)");
  45.     asm("    movl    r1,*16(ap)");
  46. }
  47.  
  48. extended_div(d, h, l, qp, rp)
  49. int d, h, l;
  50. int *qp, *rp;
  51. {
  52.     asm("    clrl    r0");
  53.     asm("    movl    8(ap),r1");
  54.     asm("    ashq    $-1,r0,r0");
  55.     asm("    addl2    12(ap),r0");
  56.     asm("    ediv    4(ap),r0,*16(ap),*20(ap)");
  57. }
  58.  
  59. #endif
  60.  
  61. #ifdef MC68K
  62.  
  63. extended_mul(d, q, r, hp, lp)
  64. int d, q, r;
  65. int *hp, *lp;
  66. {
  67.     asm("    movl    d2,a7@-");
  68.     asm("    movw    a6@(8),d0");
  69.     asm("    mulu    a6@(14),d0");
  70.     asm("    movw    a6@(10),d1");
  71.     asm("    mulu    a6@(12),d1");
  72.     asm("    addl    d1,d0");
  73.     asm("    movw    a6@(8),d1");
  74.     asm("    mulu    a6@(12),d1");
  75.     asm("    movw    a6@(10),d2");
  76.     asm("    mulu    a6@(14),d2");
  77.     asm("    swap    d2");
  78.     asm("    addw    d0,d2");
  79.     asm("    swap    d2");
  80.     asm("    swap    d0");
  81.     asm("    addxw    d0,d1");
  82.     asm("    clrl    d0");
  83.     asm("    swap    d1");
  84.     asm("    addxw    d0,d1");
  85.     asm("    swap    d1");
  86.     asm("    addl    a6@(16),d2");
  87.     asm("    addxl    d0,d1");
  88.     asm("    lsll    #1,d2");
  89.     asm("    roxll    #1,d1");
  90.     asm("    lsrl    #1,d2");
  91.     asm("    movl    a6@(20),a0");
  92.     asm("    movl    d1,a0@");
  93.     asm("    movl    a6@(24),a0");
  94.     asm("    movl    d2,a0@");
  95.     asm("    movl    a7@+,d2");
  96. }
  97.  
  98. extended_div(d, h, l, qp, rp)
  99. int d, h, l;
  100. int *qp, *rp;
  101. {
  102.     asm("    moveml    #0x3000,a7@-");
  103.     asm("    moveml    a6@(8),#0x307");
  104.     asm("    lsll    #1,d2");
  105.     asm("    addql    #1,d2");
  106.     asm("    movw    #31,d3");
  107.     asm("label2:    subl    d0,d1");
  108.     asm("    bccs    label1");
  109.     asm("    addl    d0,d1");
  110.     asm("label1:    roxll    #1,d2");
  111.     asm("    roxll    #1,d1");
  112.     asm("    dbf    d3,label2");
  113.     asm("    roxrl    #1,d1");
  114.     asm("    notl    d2");
  115.     asm("    movl    d2,a0@");
  116.     asm("    movl    d1,a1@");
  117.     asm("    moveml    a7@+,#0xc");
  118. }
  119.  
  120. #endif
  121.  
  122. #ifndef NEWS
  123. #ifdef MC68020
  124.  
  125. extended_mul(d, q, r, hp, lp)
  126. int d, q, r;
  127. int *hp, *lp;
  128. {
  129.     asm("    movl    d2,a7@-");
  130.     asm("    clrl    d2");
  131.     asm("    movl    a6@(8),d0");
  132.     asm("    mulul    a6@(12),d1:d0");
  133.     asm("    addl    a6@(16),d0");
  134.     asm("    addxl    d2,d1");
  135.     asm("    lsll    #1,d0");
  136.     asm("    roxll    #1,d1");
  137.     asm("    lsrl    #1,d0");
  138.     asm("    movl    a6@(20),a0");
  139.     asm("    movl    d1,a0@");
  140.     asm("    movl    a6@(24),a0");
  141.     asm("    movl    d0,a0@");
  142. }
  143.  
  144. extended_div(d, h, l, qp, rp)
  145. int d, h, l;
  146. int *qp, *rp;
  147. {
  148.     asm("moveml    a6@(12),#0x303");
  149.     asm("lsll    #1,d1");
  150.     asm("lsrl    #1,d0");
  151.     asm("roxrl    #1,d1");
  152.     asm("divul    a6@(8),d0:d1");
  153.     asm("movl    d0,a1@");
  154.     asm("movl    d1,a0@");
  155. }
  156.  
  157. #endif
  158. #endif
  159.  
  160. #ifdef ATT3B2
  161.  
  162. extended_mul(d, q, r, hp, lp)
  163. int d, q, r;
  164. int *hp, *lp;
  165. {
  166.     register int *r8, *r7, *r6, *r5, *r4;
  167.  
  168.     asm("    movh    0(%ap),%r5");        /* R5 <= high(d) */
  169.     asm("    andw3    &0xffff,0(%ap),%r6");    /* R6 <= low(d) */
  170.     asm("    movh    4(%ap),%r7");        /* R7 <= high(q) */
  171.     asm("    andw3    &0xffff,4(%ap),%r8");    /* R8 <= low(q) */
  172.  
  173.     asm("    mulw3    %r6,%r8,%r0");        /* R0 <= low(d)*low(q) */
  174.     asm("    andw3    &0x80000000,%r0,%r1");    /* save MSB(R0) */
  175.     asm("    lrsw3    &16,%r1,%r1");
  176.     asm("    andw2    &0x7fffffff,%r0");    /* MSB(R0) <= 0 */
  177.     asm("    addw2    8(%ap),%r0");        /* R0 <= R0 + r */
  178.     asm("    andw3    &0xffff,%r0,%r4");    /* save low(R0) */
  179.  
  180.     asm("    lrsw3    &16,%r0,%r0");        /* R0 >> 16 */
  181.     asm("    addw2    %r1,%r0");        /* resume MSB(R0) */
  182.  
  183.     asm("    mulw3    %r5,%r8,%r2");        /* R0 <= high(d)*low(q) */
  184.     asm("    addw2    %r2,%r0");        /*       + R0           */
  185.     asm("    andw3    &0x80000000,%r0,%r1");    /* save MSB(R0) */
  186.     asm("    lrsw3    &15,%r1,%r1");
  187.     asm("    andw2    &0x7fffffff,%r0");    /* MSB(R0) <= 0 */
  188.  
  189.     asm("    mulw3    %r7,%r6,%r2");        /* R0 <= low(d)*high(q) */
  190.     asm("    addw2    %r2,%r0");        /*       + R0           */
  191.     asm("    andw3    &0x7fff,%r0,%r2");    /* R2 <= low(R0) */
  192.     asm("    llsw2    &16,%r2");        /* high <= low(R2)*2^16 */
  193.     asm("    orw3    %r2,%r4,*16(%ap)");    /*         + low(R4)    */
  194.  
  195.     asm("    lrsw3    &15,%r0,%r0");        /* R0 >> 15, not 16 */
  196.     asm("    addw2    %r1,%r0");        /* resume MSB(R0) */
  197.     asm("    mulw3    %r5,%r7,%r1");        /* high <= high(d)*high(q) */ 
  198.     asm("    llsw2    &1,%r1");        /*         + R0            */
  199.     asm("    addw3    %r0,%r1,*12(%ap)");
  200. }
  201.  
  202. extended_div(d, h, l, qp, rp)
  203. int d, h, l;
  204. int *qp, *rp;
  205. {
  206.     register int *r8, *r7, *r6;
  207.  
  208.     asm("    movw    0(%ap),%r2");
  209.     asm("    movw    4(%ap),%r0");
  210.     asm("    movw    8(%ap),%r1");
  211.     asm("    movw    &0,%r8");    /* quotient to go */
  212.     asm("    movw    &31,%r7");    /* loop counter */
  213.  
  214.     asm("a:    llsw2    &1,%r8");    /* R8 << 1, no need 1st time */
  215.     asm("    llsw2    &1,%r0");    /* R0 << 1 */
  216.     asm("    llsw2    &1,%r1");    /* R1 << 1 */
  217.     asm("    jge    b");        /* skip if MSB(R1) = 0 */
  218.     asm("    orw2    &1,%r0");    /* LSB(R0) <= 1 */
  219.     asm("b:    subw3    %r2,%r0,%r6");    /* R6 <= R0 - R2 */
  220.     asm("    jl    c");        /* skip if R5 < 0 */
  221.     asm("    movw    %r6,%r0");    /* R0 <= R0 - R2 */
  222.     asm("    orw2    &1,%r8");    /* LSB(R8) <= 1 */
  223.     asm("c:    subw2    &1,%r7");    /* R7 <= R7 - 1 */
  224.     asm("    jg    a");        /* repeat while R3 > 0 */
  225.  
  226.     asm("    movw    %r8,*12(%ap)");
  227.     asm("    movw    %r0,*16(%ap)");
  228. }
  229.  
  230. #endif
  231.  
  232. #ifdef NS32K
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  
  265.  
  266.  
  267. #endif
  268.  
  269. #ifdef S3000
  270.  
  271. extended_mul(d, q, r, hp, lp)
  272. int d, q, r;
  273. int *hp, *lp;
  274. {
  275.     long long int ld, lq, lr, z;
  276.         int zh, zl;
  277.  
  278.         ld = d;
  279.         lq = q;
  280.         lr = r;
  281.         z = ld*lq+lr;
  282.         zl = (z & 0x000000007fffffffLL);
  283.         zh = (z >> 31LL);
  284.         *hp = zh;
  285.         *lp = zl;
  286. }
  287.  
  288. extended_div(d, h, l, qp, rp)
  289. int d, h, l;
  290. int *qp, *rp;
  291. {
  292.     long long int lh, ld, ll;
  293.  
  294.     ld = d;
  295.     lh = h;
  296.     ll = l;
  297.     lh = (lh << 31LL);
  298.     lh = (lh | ll);
  299.     *qp = (lh/ld);
  300.     *rp = (lh%ld);
  301. }
  302.  
  303. #endif
  304.  
  305. #ifdef IBMRT
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  
  344.  
  345.  
  346.  
  347.  
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  
  419.  
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  
  429. #endif
  430.  
  431. #ifdef NEWS
  432.  
  433. extended_mul(d, q, r, hp, lp)
  434. int d, q, r;
  435. int *hp, *lp;
  436. {
  437.     asm("    move.l    d2,-(sp)");
  438.     asm("    clr.l    d2");
  439.     asm("    move.l    (8,fp),d0");
  440.     asm("    mulu.l    (12,fp),d1:d0");
  441.     asm("    add.l    (16,fp),d0");
  442.     asm("    addx.l    d2,d1");
  443.     asm("    lsl.l    #1,d0");
  444.     asm("    roxl.l    #1,d1");
  445.     asm("    lsr.l    #1,d0");
  446.     asm("    move.l    (20,fp),a0");
  447.     asm("    move.l    d1,(a0)");
  448.     asm("    move.l    (24,a6),a0");
  449.     asm("    move.l    d0,(a0)");
  450. }
  451.  
  452. extended_div(d, h, l, qp, rp)
  453. int d, h, l;
  454. int *qp, *rp;
  455. {
  456.     asm("movem.l    (12,fp),#0x303");
  457.     asm("lsl.l    #1,d1");
  458.     asm("lsr.l    #1,d0");
  459.     asm("roxr.l    #1,d1");
  460.     asm("divu.l    (8,fp),d0:d1");
  461.     asm("move.l    d0,(a1)");
  462.     asm("move.l    d1,(a0)");
  463. }
  464.  
  465. #endif
  466.